Pacotes

library(tidyverse)
library(psych)
library(plotly)
library(GGally)
library(car)


I. Modelo Linear Normal

Para a aplicação a seguir, consideraremos o modelo de regressão normal linear

\[ y_i = \beta_1 + \beta_2x_{2i} + ... + \beta_px_{pi} + \epsilon_i, \quad i = 1,...,n \]

em que os erros \(e_i\) são variáveis aleatórias independentes, normalmente distribuídas, de média zero e variância \(\sigma^2\) constante.


Previsão do Preço de Venda de Imóveis

Dados

Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 112). Dados disponíveis em: https://www.ime.usp.br/~giapaula/textoregressao.htm

Neste exemplo, vamos modelar o preço de venda de imóveis a partir de dados relativos a uma amostra de 27 imóveis. As variáveis do conjunto de dados são:

  • imposto: imposto do imóvel (em US$ 100);

  • areaT: área do terreno (em 1.000 pés quadrados);

  • areaC: área construída (em 1.000 pés quadrados);

  • idade: idade da residência (em anos);

  • preco: preço de venda do imóvel (em US$ 1.000).

Sendo assim, o nosso objetivo é encontrar o melhor ajuste linear que explique e quantifique a variável preço de venda a partir das demais.

# Obtendo os dados
dados = read.table("dados/imoveis.dat")
colnames(dados) = c("imposto", "areaT", "areaC", "idade", "preco")
attach(dados)

# Algumas observações dos dados
head(dados)
##   imposto areaT areaC idade preco
## 1  4.9176 3.472 0.998    42  25.9
## 2  5.0208 3.531 1.500    62  29.5
## 3  4.5429 2.275 1.175    40  27.9
## 4  4.5573 4.050 1.232    54  25.9
## 5  5.0597 4.455 1.121    42  29.9
## 6  3.8910 4.455 0.988    56  29.9
# Medidas descritivas
describe(dados)
##         vars  n  mean    sd median trimmed   mad   min   max range  skew
## imposto    1 27  7.24  2.88   6.09    6.84  2.15  3.89 15.42 11.53  1.40
## areaT      2 27  6.35  2.40   5.85    6.22  2.07  2.28 12.80 10.53  0.68
## areaC      3 27  1.51  0.56   1.49    1.41  0.39  0.98  3.42  2.44  2.02
## idade      4 27 36.48 14.05  40.00   36.96 14.83  3.00 62.00 59.00 -0.34
## preco      5 27 38.50 14.31  36.90   35.65  8.90 25.90 84.90 59.00  2.25
##         kurtosis   se
## imposto     1.34 0.55
## areaT       0.01 0.46
## areaC       4.08 0.11
## idade      -0.54 2.70
## preco       4.63 2.75


Análise Descritiva

Box-plots

plot_ly(type = "box") %>%
  add_trace(y = imposto, name = "imposto") %>% 
  add_trace(y = areaT, name = "área do terreno") %>% 
  add_trace(y = areaC, name = "área construída") %>%
  add_trace(y = idade, name = "idade do imóvel") %>%
  add_trace(y = preco, name = "preço de venda")

Em um breve resumo descritivo dos dados, observamos que as altas variâncias das variáveis idade e preço destacam-se das demais. Com uma idade mediana de 40 anos, temos um imóvel de apenas 3 anos de idade e, com um preço mediano de US$ 36.900, temos dois imóveis bem mais caros nos valores de US$ 82.900 e US$ 84.900.

Os gráficos de box-plot nos mostram que os dados possuem outliers e que as variáveis possuem distribuições assimétricas. Focando na variável alvo, preço de venda, os pontos destoantes são justamente os dois imóveis mais caros.


Densidades, Dispersões e Correlações

g = ggpairs(dados, aes(color = I("slategray"), fill = I("slategray")),
            lower = list(continuous = wrap("smooth",col="black")),
            diag=list(continuous = wrap("densityDiag",alpha=0.5,size=1)))
ggplotly(g)

No gráfico acima, as curvas representadas na diagonal principal deixam mais evidente a constatação de assimetria nas distribuições observada nos box-plots. A variável preço possui correlações positivas fortes com as variáveis imposto, área do terreno e área construída, mas uma correlação negativa fraca com a variável idade. Em outras palavras, quanto maiores os valores de imposto, área construída e área do terreno, maior o preço de venda do imóvel.

Podemos observar também que há correlações relevantes entre as variáveis explicativas imposto, área construída e área do terreno. Isso nos dá indícios de multicolinearidade e, assim, uma possível redundância de informação passada por essas variáveis. Então, talvez um modelo completo com todas as variáveis não seja o mais adequado.


Modelo Completo

Nesse primeiro ajuste, consideraremos o modelo completo com todas as variáveis disponíveis.

\[ preço_i = \beta_0 + \beta_1imposto_i + \beta_2areaC_i + \beta_3areaT_i + \beta_4idade_i \]

model1 = lm(preco~.,dados)
summary(model1)
## 
## Call:
## lm(formula = preco ~ ., data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9685 -2.7152  0.2663  2.1374  6.9872 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.43587    4.09233   0.595  0.55776    
## imposto      2.07823    0.55296   3.758  0.00109 ** 
## areaT        0.23238    0.50658   0.459  0.65094    
## areaC       13.97381    2.90663   4.808  8.4e-05 ***
## idade       -0.04376    0.06628  -0.660  0.51594    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.077 on 22 degrees of freedom
## Multiple R-squared:  0.9313, Adjusted R-squared:  0.9188 
## F-statistic: 74.56 on 4 and 22 DF,  p-value: 1.809e-12

Verificamos que, para o modelo completo, as variáveis imposto e área construída são as únicas significativas e obtivemos um \(R^2\) ajustado de 0.92, muito alto, indicando um possível bom ajuste aos dados amostrais. Porém, o \(R^2\) por si só não nos garante que o modelo seja o mais adequado.


Seleção de variáveis

Aplicaremos o critério de informação de Akaike (AIC) ao modelo completo, para selecionar as variáveis que deverão permanecer.

MASS::stepAIC(model1)
## Start:  AIC=80.36
## preco ~ imposto + areaT + areaC + idade
## 
##           Df Sum of Sq    RSS    AIC
## - areaT    1      3.50 369.14 78.614
## - idade    1      7.25 372.88 78.887
## <none>                 365.64 80.357
## - imposto  1    234.76 600.40 91.748
## - areaC    1    384.13 749.77 97.746
## 
## Step:  AIC=78.61
## preco ~ imposto + areaC + idade
## 
##           Df Sum of Sq    RSS    AIC
## - idade    1     11.51 380.64 77.443
## <none>                 369.14 78.614
## - imposto  1    246.10 615.24 90.407
## - areaC    1    489.34 858.47 99.402
## 
## Step:  AIC=77.44
## preco ~ imposto + areaC
## 
##           Df Sum of Sq    RSS    AIC
## <none>                 380.64 77.443
## - imposto  1    350.01 730.65 93.049
## - areaC    1    483.16 863.80 97.569
## 
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
## 
## Coefficients:
## (Intercept)      imposto        areaC  
##      0.7903       2.2971      13.9333

O AIC nos retornou as variáveis imposto e área construída como aquelas que deverão ser mantidas no modelo, o qual chamaremos de modelo parcimonioso.


Modelo Parcimonioso

model2 = lm(preco~imposto+areaC,dados)
summary(model2)
## 
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0177 -3.3169  0.0958  2.4382  7.0949 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7903     2.2794   0.347    0.732    
## imposto       2.2971     0.4890   4.698 8.96e-05 ***
## areaC        13.9333     2.5244   5.519 1.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.982 on 24 degrees of freedom
## Multiple R-squared:  0.9285, Adjusted R-squared:  0.9225 
## F-statistic: 155.8 on 2 and 24 DF,  p-value: 1.791e-14

E, assim, obtemos um modelo em que as variáveis imposto e área construída são altamente significativas.

Porém, na análise descritiva, observamos que as variáveis área construída e área do terreno possuem uma correlação considerável. Então, podemos ainda testar um segundo modelo parcimonioso considerando a variável área do terreno ao invés da área construída.

model3 = lm(preco~imposto+areaT,dados)
summary(model3)
## 
## Call:
## lm(formula = preco ~ imposto + areaT, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.621  -3.070   0.446   3.167  10.610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.1614     3.3019   0.957   0.3479    
## imposto       3.9126     0.5293   7.391 1.24e-07 ***
## areaT         1.1017     0.6347   1.736   0.0954 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.655 on 24 degrees of freedom
## Multiple R-squared:  0.8558, Adjusted R-squared:  0.8438 
## F-statistic: 71.22 on 2 and 24 DF,  p-value: 8.081e-11

E, assim, obtemos um modelo em que a variável área do terreno passa de não significativa para pouco significativa ao nível de 10% de significância. Já a variável imposto continua bastante significativa ao nível de 1% de significância.


Diagnóstico

Agora, faremos a análise de diagnóstico dos dois modelos parcimoniosos propostos, para verificar suas adequações.

Envelope

par(mfrow=c(1,2))
qqPlot(model2, pch = 16, main = "preco ~ imposto + areaC")
## [1]  6 10
qqPlot(model3, pch = 16, main = "preco ~ imposto + areaT")

## [1]  9 27


Teste de Normalidade

\(H_0\): Os resíduos têm distribuição normal.

# Shapiro-Wilk
shapiro.test(model2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model2$residuals
## W = 0.96558, p-value = 0.4903
shapiro.test(model3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model3$residuals
## W = 0.9793, p-value = 0.8462

Ao nível de 5% de significância, não rejeitamos a hipótese de que os resíduos tenham distribuição normal. Nos gráficos de envelope para o ajuste normal, os pontos estão dentro das bandas, exceto pela observação 27 no modelo com área do terreno. Porém, a presença de leves ondulações pode ser indício de variância não constante. Vamos verificar se de fato isso ocorre no gráfico de resíduos por valores ajustados.


Resíduos x Valores ajustados

residualPlot(model2, type = "rstudent", id = TRUE,
             ylab = "resíduo studentizado", xlab = "valor ajustado",
             main = "preco ~ imposto + areaC", pch = 16, quadratic = FALSE)

residualPlot(model3, type = "rstudent", id = TRUE,
             ylab = "resíduo studentizado", xlab = "valor ajustado",
             main = "preco ~ imposto + areaT", pch = 16, quadratic = FALSE)

Embora não haja violação do pressuposto de normalidade dos resíduos, o gráfico de resíduos por valores ajustados nos dá indícios de heterocedasticidade, isto é, variância não constante, como havíamos suposto a partir do gráfico de envelope. Porém, esse comportamento pode estar sendo influenciado pelas observações destoantes vistas na etapa descritiva. Além disso, podemos perceber no gráfico de resíduos por valores ajustados a presença de pontos aberrantes, isto é, muito além do intervalo [-2,2]. No caso do primeiro modelo parcimonioso, o ponto aberrante é a observação 10, enquanto que no segundo modelo parcimonioso são as observações 9 e 27. Sendo assim, a seguir, vamos verificar se de fato esses pontos estão influenciando na qualidade dos ajustes.


Pontos Influentes

Partiremos agora para a investigação de possíveis pontos influentes, que podem estar atrapalhando a qualidade do ajuste.

Medidas de influência:

  • DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo;

  • DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores ajustados do modelo;

  • COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo;

  • HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados, é a métrica de alavancagem.

inf = influence.measures(model2)
summary(inf)
## Potentially influential observations of
##   lm(formula = preco ~ imposto + areaC, data = dados) :
## 
##    dfb.1_  dfb.imps dfb.areC dffit   cov.r   cook.d  hat    
## 9  -0.28    0.00     0.19     0.35    2.18_*  0.04    0.49_*
## 10 -1.25_*  0.29     0.59     1.62_*  0.87    0.73    0.32  
## 27 -0.13   -2.04_*   1.86_*  -2.12_*  2.04_*  1.39_*  0.61_*
inf = influence.measures(model3)
summary(inf)
## Potentially influential observations of
##   lm(formula = preco ~ imposto + areaT, data = dados) :
## 
##    dfb.1_  dfb.imps dfb.areT dffit   cov.r   cook.d  hat    
## 9  -1.15_*  1.65_*  -0.46     2.00_*  0.81    1.07_*  0.37_*
## 10 -1.21_*  0.47     0.68     1.54_*  1.02    0.69    0.35_*
## 27  0.09   -2.83_*   2.36_*  -3.05_*  0.33_*  1.85_*  0.35_*

Baseando-se nas medidas de influência acima, que verificam o efeito da remoção de cada uma das observações individualmente, temos que os possíveis pontos de influência são as observações 9 e 27, para ambos os modelos parcimoniosos, e a observação 10, para o segundo. Percebemos também que a observação 27 é a mais crítica, por impactar nas estimativas dos dois parâmetros, nos valores ajustados e na variância do modelo. Também podemos constatar isso graficamente, como se segue a baixo.

# DISTANCIA DE COOK

plot(model2, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
     sub.caption = "preco ~ imposto + areaC")

plot(model3, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
     sub.caption = "preco ~ imposto + areaT")

# ALAVANCAGEM

cut = 2*3/27 # 2*p/n
plot(hatvalues(model2), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem",
     sub = "preco ~ imposto + areaC")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model2) * as.integer(hatvalues(model2) > cut),
     cex=0.8, p=1, offset=0.3)

cut = 2*3/27 # 2*p/n
plot(hatvalues(model3), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem",
     sub = "preco ~ imposto + areaT")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model3) * as.integer(hatvalues(model3) > cut),
     cex=0.8, p=1, offset=0.3)

As observações 9, 10 e 27 são referentes aos dados:

dados[c(9,10,27),c("imposto","areaC","areaT","preco")]
##    imposto areaC areaT preco
## 9  15.4202  3.42   9.8  84.9
## 10 14.4598  3.00  12.8  82.9
## 27 12.0000  1.20   5.0  41.0

Os imóveis 9 e 10 possuem os valores de imposto mais elevados e as maiores áreas de terreno e de construção. De acordo com a relação linear identificada na etapa descritiva e a significância dessas variáveis para o modelo preditivo, podemos dizer que essas características justificam um alto preço de venda. Porém, os valores atribuídos aos imóveis 9 e 10 são muito mais elevados que os esperados pelo modelo proposto a imóveis com os mesmos atributos. Já a observação 27 causa estranheza por ter um valor de imposto tão alto e próximo aos dos imóveis 9 e 10, mesmo tendo menos da metade de suas áreas de terreno e de construção e custando pouco menos da metade que o imóvel 10.


Impacto das remoções das observações

# Compara os coeficientes e seus respectivos p-valores de um modelo
# após a remoção de uma observação
impacto = function(m1, m2){
  imp = (coef(m1) - coef(m2))/coef(m1) * 100
  impacto = data.frame(coef_com_i = coef(m1),
                       coef_sem_i = coef(m2),
                       impacto = imp,
                       pval_com_i = summary(m1)$coefficients[,4],
                       pval_sem_i = summary(m2)$coefficients[,4])
  return(impacto)
}
  • Observação i=9:

    model2_sem9 = lm(preco~imposto+areaC,dados,subset=-9)
    impacto(model2,model2_sem9)
    ##             coef_com_i coef_sem_i      impacto   pval_com_i   pval_sem_i
    ## (Intercept)  0.7903027   1.433603 -81.39918902 7.318305e-01 0.6304908219
    ## imposto      2.2970580   2.297773  -0.03113799 8.955810e-05 0.0001221436
    ## areaC       13.9332510  13.454944   3.43284823 1.122859e-05 0.0001144541
    model3_sem9 = lm(preco~imposto+areaT,dados,subset=-9)
    impacto(model3,model3_sem9)
    ##             coef_com_i coef_sem_i    impacto   pval_com_i   pval_sem_i
    ## (Intercept)   3.161412   6.558508 -107.45502 3.478823e-01 5.383568e-02
    ## imposto       3.912561   3.131158   19.97165 1.243498e-07 1.077092e-05
    ## areaT         1.101699   1.360792  -23.51762 9.543188e-02 2.722222e-02

    A remoção da observação 9 não teve impactos relevantes nos coeficientes dos dois modelos parcimoniosos e também não causou mudança inferencial, uma vez que as variáveis imposto, área construída e área do terreno permaneceram igualmente significativas nos respectivos modelos;

  • Observação i=10:

    model2_sem10 = lm(preco~imposto+areaC,dados,subset=-10)
    impacto(model2,model2_sem10)
    ##             coef_com_i coef_sem_i     impacto   pval_com_i   pval_sem_i
    ## (Intercept)  0.7903027   3.398847 -330.069021 7.318305e-01 1.640664e-01
    ## imposto      2.2970580   2.167849    5.624971 8.955810e-05 7.657851e-05
    ## areaC       13.9332510  12.571459    9.773685 1.122859e-05 2.390378e-05
    model3_sem10 = lm(preco~imposto+areaT,dados,subset=-10)
    impacto(model3,model3_sem10)
    ##             coef_com_i coef_sem_i     impacto   pval_com_i   pval_sem_i
    ## (Intercept)   3.161412  6.8915170 -117.988556 3.478823e-01 6.463495e-02
    ## imposto       3.912561  3.6816714    5.901246 1.243498e-07 2.123188e-07
    ## areaT         1.101699  0.6967438   36.757341 9.543188e-02 2.749152e-01

    A remoção da observação 10 não causou impactos relevantes nos coeficientes do primeiro modelo parcimonioso e também não ocasionou mudança inferencial, tendo as variáveis imposto e área construída mantido seus níveis de significância. Já no segundo modelo parcimonioso, a remoção da observação 10 ocasionou uma mudança inferencial, uma vez que a variável área de terreno perde sua significância;

  • Observação i=27:

    model2_sem27 = lm(preco~imposto+areaC,dados,subset=-27)
    impacto(model2,model2_sem27)
    ##             coef_com_i coef_sem_i   impacto   pval_com_i   pval_sem_i
    ## (Intercept)  0.7903027   1.068075 -35.14754 7.318305e-01 0.6320276693
    ## imposto      2.2970580   3.255657 -41.73159 8.955810e-05 0.0001915344
    ## areaC       13.9332510   9.412139  32.44836 1.122859e-05 0.0155595300
    model3_sem27 = lm(preco~imposto+areaT,dados,subset=-27)
    impacto(model3,model3_sem27)
    ##             coef_com_i coef_sem_i    impacto   pval_com_i   pval_sem_i
    ## (Intercept)   3.161412  2.9423870   6.928084 3.478823e-01 2.603227e-01
    ## imposto       3.912561  5.0694489 -29.568550 1.243498e-07 4.765596e-10
    ## areaT         1.101699 -0.0528502 104.797154 9.543188e-02 9.260592e-01

    A remoção da observação 27 impactou em mudança inferencial nos dois modelos parcimoniosos. Enquanto que no primeiro, a significância da variável área construída é reduzida, no segundo, a variável área do terreno perde completamente sua significância.


Como o segundo modelo parcimonioso não nos trouxe nenhuma melhoria para o ajuste e as conclusões sobre os pontos de influência não favoreceram a permanência da variável área do terreno, optaremos pelo primeiro modelo parcimonioso com as variáveis imposto e área construída sugerido inicialmente pelo critério de seleção de Akaike.

A partir dos gráficos de resíduos por valores ajustados, de distância de Cook e de alavancagem, podemos fazer as seguintes classificações:

  • Ponto aberrante: observação 10;

  • Pontos de alavanca: observações 9, 10 e 27;

  • Ponto influente: observação 27.


Interpretação

Então, concluímos que

\[ preço_i = 0.79 + 2.30*imposto_i + 13.93*areaC_i \quad , \]

onde:

  • A cada US$ 100,00 acrescidos no valor de imposto, aumenta-se, em média, US$ 2.300,00 no preço de venda;

  • A cada 1.000 metros quadrados de área construída acrescidos, aumenta-se, em média, US$ 13.930,00 no preço de venda.